#include "MonteCarloPi.h"

struct MonteCarloPi
{
	long long inside;	// Points generated inside
	long long total;		// Total Points generated
	long long left;		// Points still to be generated
	jkiss jk;						// random generator struct
};


void MonteCarloPi_Generate(struct MonteCarloPi *points)
{
	// Generate
	for(; points->left > 0; points->left--)
	{
		// Randomly generate points
		double x, y, dist;

		x = next(&(points->jk) );
		y = next(&(points->jk) );

		// Calculate distance to center
		dist = sqrt(x*x +y*y);

		// Increment inside and total value
		if(dist < 1)
			points->inside++;
		
		points->total++;
	
		#ifdef TESTING
		printf("\t Point:(%lf, %lf) - Distance: %lf \t [%ld Points Left | \t Thread ID: %d]\n", x, y, dist, points->left, pthread_self());
		#endif
	}

}

void MonteCarloPi_Calculate(struct MonteCarloPi values[], int size, double *pi)
{
	int i;
	long long inside = 0;
	long long total = 0;

	// Sum all up
	for(i=0; i<size; i++){
		inside += values[i].inside;
		total += values[i].total;
	}

	#ifdef TESTING
	printf("\n\t%ld - Inside Circle: ", (long)inside);
	printf("\n\t%ld - Todal        : ", (long)total);
	printf("\n\n");
	#endif

	// Calc pi = 4 *inside/total
	*pi = 4* (double) inside / (double) total;
}

double MonteCarloPi(long long iterations, int thread_count)
{
	int i;
	struct MonteCarloPi *points;	// Saved point positions
	pthread_t *thread;				// Threads
	void **status;					// Thread Status
	double pi;						// Pi result

	// Allocate
	points = (struct MonteCarloPi *) malloc(sizeof(struct MonteCarloPi) *thread_count);
	thread = (pthread_t *) malloc(sizeof(pthread_t) *thread_count);
	status = (void **) malloc(sizeof(void *) *thread_count);

	// Generate points
	for(i=0; i<thread_count; i++)
	{
		// Init random generator
		JK_Init( &(points[i].jk) );
		points[i].inside = 0;
		points[i].total = 0;
		points[i].left = (long long) (iterations/thread_count);

		// Create thread
		pthread_create(thread +i, NULL, (void * (*)(void *)) MonteCarloPi_Generate, (void *)(points +i));
	}

	// Join all threads
	for(i=0; i<thread_count; i++)
		pthread_join(thread[i], status +i);

	// Sum all and calculate pi
	MonteCarloPi_Calculate(points, thread_count, &pi);

	// Free
	free(points);
	free(thread);
	free(status);

	return pi;
}
